Hyperviscosity, Galerkin truncation and bottlenecks in turbulence 
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It is shown that the use of a high power a of the Laplacian in the dissipative term of hydrodynamical equations 
leads asymptotically to truncated inviscid conservative dynamics with a finite range of spatial Fourier modes. 
Those at large wavenumbers thermalize, whereas modes at small wavenumbers obey ordinary viscous dynamics 
[C. Cichowlas et al. Phys. Rev. Lett. 95, 264502 (2005)]. The energy bottleneck observed for finite a may be 
interpreted as incomplete thermalization. Artifacts arising from models with a > 1 are discussed. 
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A single Maxwell daemon embedded in a turbulent flow 
would hardly notice that the fluid is not exactly in thermal 
equilibrium because incompressible turbulence, even at very 
high Reynolds numbers, constitutes a tiny perturbation on 
thermal molecular motion. Dissipation in real fluids is just 
the transfer of macroscopically organized (hydrodynamic) en- 
ergy to molecular thermal energy. Artificial microscopic sys- 
tems can act just like the real one as far as the emergence 
of hydrodynamics is concerned; for instance, in lattice gases 
the "molecules" are discrete Boolean entities |Q]] and ther- 
malization is easily observed at high wavenumbers J2l. An- 
other example has been found recently by Cichowlas et al. ^ 
wherein the Euler equations of ideal non-dissipative flow are 
(Galerkin) truncated by keeping only a finite - but large - 
number of spatial Fourier harmonics. The modes with the 
highest wavenumbers k then rapidly thermalize through a 
mechanism discovered by T.D. Lee 0] and studied further 
by R.H. Kraichnan [5], leading in three dimensions (3D) to 
an equipartition energy spectrum oc k 2 . The thermalized 
modes act as a fictitious microworld on modes with smaller 
wavenumbers in such a way that the usual dissipative Navier- 
Stokes dynamics is recovered at large scales [25]. 

All the known systems presenting thermalization are con- 
servative. As we shall show themalization may be present in 
dissipative hydrodynamic systems when the dissipation rate 
increases so fast with the wavenumber that it mimics ideal 
hydrodynamics with a Galerkin truncation. This is best un- 
derstood by considering hydrodynamics with hyperviscosity: 
the usual momentum diffusion operator (a Laplacian) is re- 
placed by the ath power of the Laplacian, where a > 1 is the 
dissipativity. Hyperviscosity is frequently used in turbulence 
modeling to avoid wasting numerical resolution by reducing 
the range of scales over which dissipation is effective [6]. 

The unforced hyperviscous ID Burgers and multi- 
dimensional incompressible Navier-Stokes (NS) equations 
are: 

d t v + vd x v = -pk G - 2a {-d 2 x ) a v- (i) 

dtV + v-Vv = -Vp-pk c r 2a (-V 2 ) a v; V-t> = . (2) 



The equations must be supplemented with suitable initial and 
boundary conditions. We employ 27r-periodic boundary con- 
ditions in space, so that we can use Fourier decompositions 
such as v(x) = J^k G lk x - Note that minus the Laplacian 
is a positive operator, with Fourier transform fc 2 , which can 
be raised to an arbitrary power a. The coefficient /i is taken 
positive to make the hyperviscous operator dissipative. The 
Galerkin wavenumber k G > is chosen off-lattice so that no 
wavenumber is exactly equal to k G . In Fourier space the hy- 
perdissipation rate is p{k j k G ) 2a , where k = \k\. 

If we now hold /i and k G fixed and let a — > oo we see 
that the hyperdissipation rate tends to zero, for k < k G , and 
to infinity, for k > k G . This implies that in the limit of infi- 
nite dissipativity, the solution of a hyperviscous hydrodynam- 
ical equation converges to that of the corresponding inviscid 
equations Galerkin-truncated at wavenumber k G . 

To define inviscid Galerkin truncation precisely, we rewrite 
Eqs. (Q]l and (0 in the abstract form d t v — B(v,v) + L a v, 
where B is a quadratic form representing the nonlinear term 
(including the pressure p in the NS case). The truncation pro- 
jector Pk G is the linear, low-pass filtering operator that, when 
applied to v, sets all Fourier harmonics with k > k G to zero. 
The inviscid, Galerkin-truncated equation, with initial condi- 
tion vq, is 

d t u = P kG B(u,u), u = Pk G v - (3) 

Since u can be written in terms of a finite number of modes 
with k < k G , Eq. (O is a dynamical system of finite dimen- 
sion. In addition to momentum, it conserves the energy and 
other quadratic invariants for the inviscid equations [5]. There 
is good numerical evidence - but no rigorous proof - that the 
solutions of the Galerkin-truncated inviscid Burgers and 3D 
Euler equations tend, at large times, to statistical equilibria 
defined by their respective invariants. 

A rigorous proof of the convergence, as a — > oo, of solu- 
tions of the hyperviscous Burgers equation CO and of the hy- 
perviscous NS equation (|2]) in any dimension to those of the 
associated Galerkin-truncated, inviscid equation will be given 
elsewhere. It uses standard tools of functional analysis; note 
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that the formidable mathematical difficulties that beset the or- 
dinary (a = 1) 3D NS equation disappear for a > 5/4 ||7|]. 

From a physicist's point of view the convergence result 
looks rather obvious, though it has hardly been noted before 
(see, however, Refs. H u\idj\): as a — > oo all the modes 
with k > k G are immediately suppressed by an infinite dissi- 
pation, whereas those with k < k G obey inviscid truncated 
dynamics. Not surprisingly, the fate of couplings between 
triads of modes whose wavenumbers straddle k G is a deli- 
cate point. In a Galerkin truncation any such triad should 
be left out. It may be shown that for a — > oo such strad- 
dling couplings are suppressed, not only for the Burgers and 
NS equation but also for the hyperviscous magnetohydrody- 
namical equations and for some turbulence closures, specifi- 
cally, the Direct Interaction Approximation (DIA) [8] and the 
Eddy-Damped-Quasi-Normal-Markovian (EDQNM) approx- 
imation 111 ill . Hence the convergence to the corresponding 
Galerkin-truncated equations holds for all the aforementioned 
equations in any dimension of space. 

There are, however, interesting exceptions among hydrody- 
namical equations for which the result does not hold. They in- 
clude the kinetic theory of resonant wave interactions 111 211 and 
the Markovian Random Coupling Model [13]. Indeed, the res- 
onant wave interaction theory arises in the limit when the pe- 
riod of the waves goes to zero and this limit does not commute 
with the limit of a vanishing damping time for modes having 
k > k G ; a similar remark can be made about the MRCM equa- 
tion. 



Let us stress that systems with a finite dissipativity - how- 
ever large - are quite different from Galerkin-truncated sys- 
tems. For example, consider the 3D NS equation with a ran- 
dom force, delta-correlated in time, for which we know the 
mean energy input e per unit volume. It is still true that, for 
a — ► oo, the solution of this equation converges to that of 
the Galerkin-truncated equation, but this time with a random 
force. If Eq is the initial energy, this solution has a mean en- 
ergy E(t) — Eq + et, which grows indefinitely in time. But, 
as soon as a is given a finite value, however large, a statistical 
steady state, in which energy input and hyperviscous energy 
dissipation balance, is achieved at large times. Such a steady 
state presents an interesting interplay of thermalization and 
dissipation, when a is large, as we show below. 

The direct numerical simulation (DNS) of the Galerkin- 
truncated 3D Euler equations in Ref. used 1600 3 Fourier 
modes. Large-a simulations of Eq. d2|i would require signifi- 
cantly higher resolution to identify the various spectral ranges 
that we can expect, namely, inertial, thermalized, and far- 
dissipation ranges and transition regimes between these. For- 
tunately, Bos and Bertoglio 11411 have shown that key features 
of the Galerkin-truncated Euler equations, such as the pres- 
ence of inertial and thermalized ranges, can be reproduced 
by the two-point EDQNM closure [11] for the energy spec- 
trum. For Eq. (O, with stochastic, white-in-time, homoge- 
neous and isotropic forcing with spectrum F(k), the hyper- 
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Here E(k,t) is the energy spectrum, At defines the set of 
p > and q > such that k, p, q can form a triangle, x, 
y, z are the cosines of its angles and the eddy-damping pa- 
rameter A is expressed in terms of the Kolmogorov constant. 
The EDQNM equations have been studied numerically for 
more than three decades lfl5ll . but their hyperviscous versions 
Eq. (HI present new difficulties that we overcome as follows. 
Since we are interested in the steady state we use an iterative 
method: the emission term, E(p, t)E(q, t) in Eq. (O, is con- 
sidered as a renormalization of the force F(k); the absorption 
term, E(q, t)E(k, t), is treated as a renormalization of the hy- 
perviscous damping [16]. We then construct a sequence of 
energy spectra that, at stage n + 1, is just the renormalized 
force divided by the renormalized damping, both based on 
stage n. This gives rapid convergence to the steady state at low 
wavenumbers; but, beyond a certain (a-dependent) wavenum- 
ber, convergence slows down dramatically and it is better to 
use time marching to obtain the steady state. At large val- 
ues of k and a the problem becomes very stiff, so we use a 
slaved fourth-order Runge-Kutta scheme IU7I1 . We discretize 
k logarithmically, with N c collocation points per octave. Triad 
interactions involving wavenumber ratios significantly larger 
than N c are poorly represented 11811 : so, since wavenumber 
ratios of up to 50 play an important role for large a, we have 
used iV c < 90; this is computationally demanding because 
the complexity of the code is 0(N 3 ). We force at the low- 
est wavenumber (k = 1) in our numerical study of Eq. © 
with k G = 10 5 , A = 0.36, and 1 < a < 729. The re- 
sulting compensated, steady-state energy spectra k +5 ^ 3 E(k) 
are shown in Fig. [TJ flat regions, extending over two to five 
decades of k (depending on a), are close to the Kolmogorov 
inertial range; for large a there is a distinct thermalized range 
with E(k) ~ k 2 (also found in the transition between classical 
and quantum superfluid turbulence [19]), as we expect from 
our discussion of the Galerkin-truncated Euler equations j26ll . 
In the far-dissipation range k > k G the spectra fall off very 
rapidly M27I1 . For all values of a the far-dissipation range is 
preceded by a bump or bottleneck. It is also observed, in some 
experiments l!20ll and DNS of Navier-Stokes, with a shape 
that is quite independent of the Reynolds number 12111 . The 
bottleneck for a = 1 has previously been explained as the in- 
hibition of the energy cascade from low to high wavenumbers 
because of viscous suppression of the cascade in the dissipa- 
tion range M22I1 . Our work provides an alternative explanation: 
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use 2 14 collocation points and set /i = 1, k G = 342.1 and 
a = 1000. In Fig. [2] we show the Burgers energy spectrum 



FIG. 1: (Color online) Log-log plots of the compensated spectrum 
k 5 ^ 3 E(k) versus k from a numerical integration of the hyperviscous 
EDQNM Eq. (O for different values of a; inset: enlarged spectra 
showing a secondary bottleneck (see text). 



the usual bottleneck may be viewed as incomplete thermaliza- 
tion. 

At large values of a the thermalized range gives rise to an 
eddy viscosity v c ddy This acts on modes with wavenumbers 
lower than those in the thermalized range; the corresponding 
damping rate is v e ddyk 2 . The eddy viscosity can be expressed 
as an integral over the thermalized range II14I18H. As a grows, 
so does ^eddy and, eventually, the renormalized viscous damp- 
ing overwhelms the hyperviscous damping for modes at low 
wavenumbers (below those in the thermalized range). The dy- 
namics of these modes is then governed by the usual a = 1 
equation. Not surprisingly, then, we find a pseudo-dissipation 
range around k ~ 10 4 that is shown in an expanded scale in 
the inset of Fig. [2 a similar range for the Galerkin-truncated 
case is discussed in Ref. B1411 and is already visible in the DNS 
of Ref. Jit]. For large a the inset of Fig. [TJ also shows a sec- 



ondary bottleneck range for 10 < k < 10 ; this may be 
viewed as the usual (a = 1) EDQNM bottleneck stemming 
fromj/cddy H. 

Our results apply to compressible flows also. We have stud- 
ied the simplest instance, that is the unforced hyperviscous 
ID Burgers equation ([TJ [29]. Its solution converges to the 
entropy solution, i.e., the standard solution with shocks, ob- 
tained when k G — > oo for any a > 1 I23I1 . Here we are inter- 
ested in the large-a behavior at fixed k G . We do not have to 
resort to closure now since we can solve the primitive equation 
<[TJ directly by a pseudospectral method. If we choose a sin- 
gle initial condition the resulting spectrum is noisy because, 
unlike the ordinary Burgers equation, its Galerkin-truncated 
version and thus also the high-a versions are believed to be 
chaotic dynamical systems [24]. So we solve (TJ with the two- 
mode random initial condition vq(x) = s'mx + sin(2.T + <f>), 
where (f> is distributed uniformly in the interval [— tt, it]. We 
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FIG. 2: (Color online) Log-log plots of averaged energy spectra (see 
text) versus wavenumber k for the hyperviscous Burgers equation l[T) 
at various times. 

E(k) =| v(k) | 2 , averaged over 20 realizations of the phase 
4> at various times. At the latest output times the spectrum is 
almost completely flat, i.e., thermalized, with equipartition of 
the energy between all the Fourier modes. At earlier times 
E(k) behaves approximately as fc~ 2 in an inertial range that 
corresponds to shocks in physical space; there is a thermal- 
ized range at higher wavenumbers up to fc G ; for k > k G the 
spectrum falls very rapidly. No pseudo-dissipation range is 
observed here between the inertial and thermalized ranges as 
seen in the 3D NS case (Fig. [TJ. Perhaps the data are too 
noisy, but a careful examination of v(x) in physical space in- 
dicates that this phenomenon might arise from the compress- 
ible nature of the Burgers dynamics: thermalization begins 
over the whole physical range (as high-frequency noise with 
wavenumber m fc G ); noise generated close to shocks is ab- 
sorbed by them and not enough is left to produce any appre- 
ciable eddy viscosity that could broaden the shocks. 

We now summarize our main findings from the study of 
hyperviscous hydrodynamical equations with powers a of the 
Laplacian ranging from unity to very large values. 

The simplest results are obtained for very large a. The 
solutions of the ID Burgers equation or the Navier-Stokes 
equations in any space dimension d are then very close to the 
solutions of the corresponding Galerkin-truncated equations, 
displaying thermalization at wavenumbers below k G . The de- 
tailed scenario will of course be affected by the dimension 
of space. In 3D, with enough resolution, we may be able to 
observe up to five ranges: an inertial range, a secondary bot- 
tleneck, a pseudo-dissipation range, a thermalized range, and 
a far dissipation range. Because of enstrophy conservation 
and of the predominance of Fourier-space nonlocal interac- 
tions, the 2D case is rather special and deserves a separate 
study & 
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The most relevant case is of course that of ordinary dissi- 
pation (a — 1). The energy-spectrum bottleneck generally 
observed at high Reynolds numbers in 3D incompressible tur- 
bulence may be viewed as an incomplete thermalization: as 
we increase a larger and larger bottlenecks are present, even- 
tually displaying thermalization on their rising side. 

We finally deal with the case of moderately large a of the 
sort used in many simulations [6]. How safe is this procedure 
and what kind of artifacts can we expect? 

Using large values of a in simulations to "avoid wasting 
resolution" is hardly advocated by anybody, but we now un- 
derstand what goes wrong: a huge thermalized bottleneck 
will develop at high wavenumbers, whose action on smaller 
wavenumbers is an ordinary a = 1 dissipation with an eddy 
viscosity much larger than what would be permissible in a 
normal a = 1 simulation. 

When a is chosen just a bit larger than unity (e.g. a — 
2 which is standard in oceanography ||6|]) the advantage of 
widening the inertial range may be offset by artifacts at bot- 
tleneck scales; indeed, even an incomplete thermalization will 
bring the statistical properties of such scales closer to Gaus- 
sian, thereby reducing the rather strong intermittency which 
would otherwise be expected [31]. For similar reasons spu- 
rious isotropization can be expected for problems with an 
anisotropic constraint, such as rapidly rotating or stratified 
flow or MHD with a strong uniform magnetic field. 
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